function [f0,Q,fmin,fmax] = resonance_peak_properties(f,G)
  
  % given amplitude spectrum G(f), find maximum amplitude resonance peak and
  % return frequency f0 and quality factor Q of that peak
  
  % classic method: find max and full width at half max
  
  Gmax = max(G);
  
  i = find(G==Gmax); % index of max G
  % (this might need to be modified if there are several resonance peaks--
  % probably you would want to only pass some part of the spectrum to this code)
  
  f0 = f(i); % resonance frequency
  
  % find indices of frequencies at half max amplitude
  % (this part will require modification if
  % noise level is higher than half max amplitude)
  imin = find(G(1:i  )<=0.5*Gmax,1,'last');
  imax = find(G(i:end)<=0.5*Gmax,1,'first')+i-1;
  
  fmin = f(imin);
  fmax = f(imax);
  
  Q = f0/(fmax-fmin); % quality factor
  
  % alternative method using least-squares fitting to expansion around resonance
  % M. P. Robinson and J. Clegg, IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, VOL. 47, NO. 2, MAY 2005
  % this method should be more robust, so is recommended over classic method above
  
  %plot(f(imin:imax),1./G(imin:imax)),pause % uncomment to plot
  % the shape should look like parabola, and the method is to do a linear least-squares fit
  % of a parabola to the data plotted in this manner
  
  p = polyfit(f(imin:imax),1./G(imin:imax),2); % parabolic fit
  
  % extract resonance frequency and quality factor from fit,
  % using coefficients p of parabola
  f0_fit = -p(2)/(2*p(1));
  Q_fit = 0.5/sqrt(4*p(1)*p(3)/p(2)^2-1);
  
  % uncomment to compare two methods
  %disp(['resonance frequency: ' num2str([f0 f0_fit])])
  %disp(['quality factor:      ' num2str([Q Q_fit])])
  
  f0 = f0_fit; Q = Q_fit; % return values from least-squares fitting method